#Replication File for "Strategic Ambiguity, Strategic Clarity, and Dual Clarity"
#Wang et al. 2024 on Foreign Policy Analysis

rm(list=ls(all=TRUE))

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


library(foreign)
library(car)

#Please put the sav file in the same folder of this r code file

fullpath=rstudioapi::getSourceEditorContext()$path   #Make sure where is your working file location
currentfilepath=substr(fullpath,1,nchar(fullpath)-40) #49=length of r code name (including.r)
currentfilepath
setwd(currentfilepath)

#Import dataset from PP2197EC2.sav
project1 = read.spss("PP2197EC2.sav",to.data.frame=TRUE,use.value.labels=FALSE)

#Respondent's attitude on independence before the treatments
project1$ind_pre=recode(project1$Q16,"1=2;2=1;3=-1;4=-2;else=0")

#P.10, pre-attitude on independence, 10.2% oppose strongly, 19.5% support very much
prop.table(table(project1$ind_pre))

#Experimental Assignment

project1$group3=NA
project1$group3[project1$Q17a<91]="US promise"     
project1$group3[project1$Q17b<91]="US not sure"    
project1$group3[project1$Q17c<91]="US limits ind"  
table(project1$group3)

#Respondent's attitude on independence after the treatments

project1$ind_post=NA
project1$temp=NA

project1$temp[project1$Q17a<91]=project1$Q17a[project1$Q17a<91]
project1$temp[project1$Q17b<91]=project1$Q17b[project1$Q17b<91]
project1$temp[project1$Q17c<91]=project1$Q17c[project1$Q17c<91]
table(project1$temp)
project1$ind_post=recode(project1$temp,"1=2;2=1;3=-1;4=-2;else=0")


#=====Randomization Test=====#

summary(aov(ind_pre~group3,data=project1))  
#p=0.35, p.12, no sig difference on supporting independence before the treatment

#ANOVA for control variables, and Appendix Table A2

#Gender
table(project1$sex)
project1$female=NA
project1$female=recode(project1$sex,"1=0;2=1")

#Age
table(project1$age)

#Level of Education
table(project1$edu)

#Taiwan Identity 
table(project1$Q26)
project1$TWID=recode(project1$Q26,"1=1;else=0")

#Partisan Identity
table(project1$Q27)
project1$BLUE=recode(project1$Q27,"1=1;3=1;5=1;6=1;else=0")
project1$GREEN=recode(project1$Q27,"2=1;4=1;7=1;8=1;9=1;else=0")

project1$NONE=1
project1$NONE[project1$BLUE==1]=0
project1$NONE[project1$GREEN==1]=0
table(project1$NONE)
summary(project1)

#ANOVA Test on covariates (p.12)
summary(aov(sex~group3,data=project1))  #p = 0.03, sig
summary(aov(age~group3,data=project1))  #p < 0.01, sig  
summary(aov(edu~group3,data=project1))  #p = 0.76, not sig
summary(aov(TWID~group3,data=project1)) #p = 0.27, not sig
summary(aov(BLUE~group3,data=project1)) #p = 0.12, not sig
summary(aov(GREEN~group3,data=project1))#p = 0.52, not sig

#=====Results=====#

#Pre-treatment attitude (p.12)
summarySE(data=project1,"ind_pre","group3")
#DC: 0.298
#SA: 0.454
#SC: 0.375
summary(aov(ind_pre~group3,data=project1))   #ANOVA not sig


#Post-treatment attitude (p.12)
summarySE(data=project1,"ind_post","group3")
#DC: 0.056
#SA: 0.171
#SC: 0.485
summary(aov(ind_post~group3,data=project1))  #ANOVA p=0.000333 <0.01


#Figure 1 (p.13)
plot1=summarySE(data=project1,"ind_pre","group3")
plot2=summarySE(data=project1,"ind_post","group3")
plot3=NULL
plot3$group3=rep(plot1$group3,2)
plot3$ind=c(plot1$ind_pre,plot2$ind_post)
plot3$ci=c(plot1$ci,plot2$ci)
plot3$gp=c(rep("Pre-treatment",3),rep("Post-treatment",3))

plot3=data.frame(plot3)
plot3

plot3$gp=relevel(as.factor(plot3$gp), ref="Pre-treatment")
library(ggplot2)
ggplot(data=plot3,aes(x=gp,colour=group3,shape=group3,y=ind,group=group3))+
  geom_errorbar(aes(ymin=ind-ci, ymax=ind+ci),  width=.2, position=position_dodge(0.4)) +
  geom_line(position=position_dodge(0.4)) +
  geom_point(position=position_dodge(0.4), size=3)+
  xlab("")+ylab("Level of independence support (-2 to +2)")+theme_bw()+
  theme(legend.position = "none") +ylim(c(-0.1,0.8))+
  geom_text(x=2.16,y=0.5,label="US Promise (SC)\nPaired T-test p < 0.001\n(n=301)",colour="blue",hjust = 0,size=2.5)+
  geom_text(x=2.035,y=0.18,label="US Not Sure (SA)\nPaired T-test p < 0.001\n(n=304)",colour="green4",hjust = 0,size=2.5)+
  geom_text(x=1.9,y=-0.05,label="Dual Clarity (DC)\nPaired T-test p < 0.001\n(n=305)",colour="red",hjust = 0,size=2.5)+
  geom_segment(aes(x=0.8, xend=1.2, y=0.73, yend=0.73),colour="black") + 
  geom_segment(aes(x=1.8, xend=2.2, y=0.73, yend=0.73),colour="black") + 
  geom_text(x=1,y=0.755,label="ANOVA, p = 0.35",colour="black",size=2.5)+
  geom_text(x=2,y=0.755,label="ANOVA, p < 0.001",colour="black",size=2.5)


#Paired t-test (p.13)

t.test(project1$ind_pre[project1$group3=="US promise"],project1$ind_post[project1$group3=="US promise"],paired=TRUE)   #SC
t.test(project1$ind_pre[project1$group3=="US not sure"],project1$ind_post[project1$group3=="US not sure"],paired=TRUE)  #SA
t.test(project1$ind_pre[project1$group3=="US limits ind"],project1$ind_post[project1$group3=="US limits ind"],paired=TRUE)  #DC


#Robustness: Binary Coding (p.13, Appendix Figure A4)


project1$ind_post2_2=recode(project1$temp,"1=1;2=1;3=0;4=0;else=0")
table(project1$ind_post2_2,project1$group3)
prop.table(table(project1$ind_post2_2,project1$group3),2)

project1$ind_pre2_2=recode(project1$Q16,"1=1;2=1;3=0;4=0;else=0")
table(project1$ind_pre2_2,project1$group3)
prop.table(table(project1$ind_pre2_2,project1$group3),2)



project1$diff_2=project1$ind_post2_2-project1$ind_pre2_2

summarySE(data=project1,"diff_2","group3")

summary(aov(diff_2~group3,data=project1))

t.test(project1$ind_pre2_2[project1$group3=="US promise"],project1$ind_post2_2[project1$group3=="US promise"],paired=TRUE)
t.test(project1$ind_pre2_2[project1$group3=="US not sure"],project1$ind_post2_2[project1$group3=="US not sure"],paired=TRUE)
t.test(project1$ind_pre2_2[project1$group3=="US limits ind"],project1$ind_post2_2[project1$group3=="US limits ind"],paired=TRUE)


plot1=summarySE(data=project1,"ind_pre2_2","group3")
plot2=summarySE(data=project1,"ind_post2_2","group3")
plot3=NULL
plot3$group3=rep(plot1$group3,2)
plot3$ind=c(plot1$ind_pre2_2,plot2$ind_post2_2)
plot3$ci=c(plot1$ci,plot2$ci)
plot3$gp=c(rep("Pre-treatment",3),rep("Post-treatment",3))

plot3=data.frame(plot3)
plot3

plot3$gp=relevel(as.factor(plot3$gp), ref="Pre-treatment")
library(ggplot2)
ggplot(data=plot3,aes(x=gp,colour=group3,shape=group3,y=ind,group=group3))+
  geom_errorbar(aes(ymin=ind-ci, ymax=ind+ci),  width=.2, position=position_dodge(0.4)) +
  geom_line(position=position_dodge(0.4)) +
  geom_point(position=position_dodge(0.4), size=3)+
  xlab("")+ylab("Level of independence support (0 or 1)")+theme_bw()+
  theme(legend.position = "none") +ylim(c(0.4,0.8))+
  geom_text(x=1.9,y=0.6641,label="US Promise (SC)\nPaired T-test p =0.406\n(n=301)",colour="red",hjust = 0,size=2.5)+
  geom_text(x=2.17,y=0.572,label="US Not Sure (SA)\nPaired T-test\n p < 0.001(n=304)",colour="blue",hjust = 0,size=2.5)+
  geom_text(x=1.9,y=0.46,label="Dual Clarity (DC)\nPaired T-test p < 0.001\n(n=305)",colour="green4",hjust = 0,size=2.5)+
  geom_segment(aes(x=0.8, xend=1.2, y=0.73, yend=0.73),colour="black") + 
  geom_segment(aes(x=1.8, xend=2.2, y=0.73, yend=0.73),colour="black") + 
  geom_text(x=1,y=0.755,label="ANOVA, p = 0.35",colour="black",size=2.5)+
  geom_text(x=2,y=0.755,label="ANOVA, p < 0.001",colour="black",size=2.5)


#==========Table 1================
project1$group3=relevel(as.factor(project1$group3), ref="US promise")

model1=lm(ind_post~group3,data=project1)
model2=lm(ind_post~group3+ind_pre,data=project1)
model3=lm(ind_post~group3+ind_pre+female+age+edu,data=project1)
model4=lm(ind_post~group3+ind_pre+female+age+edu+TWID+BLUE+GREEN,data=project1)
library(stargazer)

stargazer(model1,model2,model3,model4,type="text")

#Robustness Check, Ordinal Regression, Appendix Table A3
library(MASS)
model1=polr(as.factor(ind_post)~group3,data=project1)
model2=polr(as.factor(ind_post)~group3+ind_pre,data=project1)
model3=polr(as.factor(ind_post)~group3+ind_pre+female+age+edu,data=project1)
model4=polr(as.factor(ind_post)~group3+ind_pre+female+age+edu+TWID+BLUE+GREEN,data=project1)

stargazer(model1,model2,model3,model4,type="text")

#Robustness Check, Table 2, p.17

project1$group3_r=NA
project1$group3_r[project1$group3=="US promise"]="SC"
project1$group3_r[project1$group3=="US limits ind"]="DC"
project1$group3_r[project1$group3=="US not sure"]="SA"
project1$group3_r=relevel(as.factor(project1$group3_r), ref="SC")

#Preceived Threat from China
project1$war=recode(project1$Q14,"1=1;2=0;3=-1")
table(project1$Q11)
#Belief in the United States
project1$US=project1$Q11
project1$US[project1$US>10]=5


model5=lm(ind_post~group3_r*war+TWID+ind_pre+female+age+edu+TWID+GREEN+BLUE,data=project1)
summary(model5)
model6=lm(ind_post~group3_r*US+TWID+ind_pre+female+age+edu+TWID+GREEN+BLUE,data=project1)
summary(model6)
stargazer(model5,model6,type="text")

#Calculating the Robust standard error, Appendix Table A5
library(lmtest)
library(sandwich)

coeftest(model5, vcov = vcovHC(model5, type = 'HC0'))
coeftest(model6, vcov = vcovHC(model6, type = 'HC0'))


#==========Figure 2=================#
library(sjPlot)

plot_model(model5, type = "pred", terms = c("war","group3_r"))+theme_bw()+
  xlab("Perceived military threat from CN (-1 to 1)")+
  ylab("Simulated level of post-treatment ind support")+ggtitle("")+
  theme(legend.title = element_blank())

#==========Figure 3=================#

plot_model(model6, type = "pred", terms = c("US","group3_r"))+theme_bw()+
  xlab("Level of Trust between CN and US (CN 0 to US 10)")+
  ylab("Simulated level of post-treatment ind support")+ggtitle("")+
  theme(legend.title = element_blank())



